library(tidyverse)
library(terra)
library(sf)
library(leaflet)Data for testing IHS
Obtain DEM
Pulling out a patch of publicly available LiDAR near Levin, Manawatū-Whanganui, on New Zealand’s North Island, from the LINZ data portal.
This dataset can help me identify the tiles I want. This has been downloaded as a geopackage and unzipped to data_spatial/aoi.
tiles <-
read_sf(file.path('data_spatial', 'aoi',
'manawatu-whanganui-lidar-1m-index-tiles-2015-2016.gpkg'))
# had a gander in QGIS, decided I want this bit:
wishlist <- tiles %>%
dplyr::filter(TileName %in%
c('BN33_1000_1444', 'BN33_1000_1445',
'BN33_1000_1544', 'BN33_1000_1545'))
aoi <- st_as_sf(st_union(wishlist))The data itself is downloadable in GeoTiff format from this dataset.
Unfortunately the dataset only has view services (WMTS and XYZ) and a point-data query API, so the tiles have to be downloaded manually from the ‘tiles table’ interface, e.g.:
This has been done and the four tiles unzipped and placed in data_spatial\elevation
elev_dir <- file.path('data_spatial', 'elevation')
list.files(elev_dir, pattern = 'DEM.*\\.tif$')
#> [1] "DEM_BN33_2015_1000_1444.tif" "DEM_BN33_2015_1000_1445.tif"
#> [3] "DEM_BN33_2015_1000_1544.tif" "DEM_BN33_2015_1000_1545.tif"The easiest way to start working with the tiles is to create a VRT over them using GDAL. Below, my OSGeo4W GDAL installation is accessed to do this.
if(!file.exists(file.path('helpers', 'tile_files.txt'))) {
tile_files <-
list.files(elev_dir, pattern = 'DEM.*\\.tif$', full.names = TRUE)
txt_con <- file(file.path('helpers', 'tile_files.txt'))
writeLines(tile_files, txt_con)
close(txt_con)
}
source(file.path('helpers', 'gdal_env.R')) # environment variables
system2('gdalbuildvrt',
args =
c('-input_file_list', file.path('helpers', 'tile_files.txt'),
file.path(getwd(), elev_dir, 'lidar_1m.vrt')))This can now be opened in R just like any other raster data source:
dem_file <- file.path(getwd(), elev_dir, 'lidar_1m.vrt')
dem <- terra::rast(dem_file)
dem
#> class : SpatRaster
#> dimensions : 1440, 960, 1 (nrow, ncol, nlyr)
#> resolution : 1, 1 (x, y)
#> extent : 1800640, 1801600, 5503200, 5504640 (xmin, xmax, ymin, ymax)
#> coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (EPSG:2193)
#> source : lidar_1m.vrt
#> name : lidar_1mA hillshade would also come in handy,
system2('gdaldem',
args =
c('hillshade', '-multidirectional',
dem_file,
file.path(getwd(), elev_dir, 'hillshade_1m.tif'))
)
hs <- rast(file.path(getwd(), elev_dir, 'hillshade_1m.tif'))